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Abstract 


MONTE CARLO MODELS AND ANALYSIS OF 
GALACTIC DISK GAMMA-RAY BURST DISTRIBUTIONS 


Jon Hakkila 

Assistant Professor of Astronomy 
Department of Mathematics/ Astronomy/ and Statistics 
Mankato State University 
MankatO/ MN 


Gamma-ray bursts are transient astronomical phenomena which have no 
quiescent counterparts in any region of the electromagnetic spectrum. 
Although temporal and spectral properties indicate that these events 
are likely energetic/ their unknown spatial distribution complicates 
astrophysical interpretation . 

Monte carlo samples of gamma-ray burst sources are created which 
belong to Galactic disk populations. Spatial analysis techniques are 
used to compare these samples to the observed distribution. From 
this, both quantitative and qualitative conclusions are drawn 
concerning allowed luminosity and spatial distributions of the actual 
sample . 

Although the BATSE experiment on GRO will significantly improve 
knowledge of the gamma-ray burst source spatial characteristics within 
only a few months of launch, the analysis techniques described herein 
will not be superceded. Rather, they may be used with BATSE results 
to obtain detailed information about both the luminosity and spatial 
distributions of the sources. 
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Introduction 


A. General Properties of Gamma-Ray Burst Sources 


Since their discovery via examination of Vela satellite data records 
(Klebesadel, Strong, and Olson 1973), gamma-ray burst sources have 
remained one of the most enigmatic classes of objects in modern 
astronomy. The strong bursts of observed gamma-radiation that these 
objects emit have no quiescent counterparts at any other wavelength of 
the electromagnetic spectrum, and are at present of unknown origin. 
Multi-wavelength observations of the transient events are also 
unverified, although sporadic visible flashes have been reported (e.g. 
Schaefer 1981) . 

The burst durations range from I < 0.1 sec to T > 100 sec, although the 
majority apparently lie between 3 and 20 seconds (Hurley, 
unpublished) . The events span a wide range of relative rise and decay 
times (Barat et al. 1984), although their general rapidity, as well as 
variations on timescales as short as milliseconds, suggest that at 
least some of the events must be associated with compact objects 
(r < 10 3 km) . Since the events are time-integrated, the registered 
output is generally measured in fluence (erg cm“^) instead of more 
conventional flux units (erg cm"^ sec“^) . 

Although the peak burst power output lies between 150 and 500 keV for 
most bursts, the distribution (as observed from the KONUS catalogue) 
is sharply peaked around 200 keV (Higdon and Lingenfelter 198 6) . SMM 
observed hard high-energy tails, showing that burst power decreases at 
high energies (Nolan et al. 1984a) . Burst power is also observed to 
turn over at low energies (e.g. Katoh et al. 1984) . There are 
indications of cyclotron absorption features in the spectra from KONUS 
(Mazets et al. 1981), HEAO-A4 (Hueter et al. 1984), and GINGA 
(Murakami et al. 1988, Fenimore et al . 1989) . Additionally, 
controversial observations of features thought to be positron-electron 
annihilation lines in emission have been made by KONUS (Mazets et al, 
1981) and by SMM (Nolan et al 1984b) . 
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B. Spatial Analyses of Burst Sources 


The location of the burst sources in space is difficult to determine. 
Without prior knowledge of their luminosities, their observed fluxes 
cannot accurately be converted into spatial positions. Astrophysical 
models for the bursts must therefore remain vague until this 
controversy is resolved, as the burst sources might be local, 
disk-population Galactic, halo-population Galactic, or extragalactic 
in origin. There are two approaches by which this spatial 
distribution can be studied; (1) by examining the angular distribution 
of sources, (2) by examining the radial distribution of sources (such 
as is attempted using the log (N) -log (S) method or the V/V max test). 


1. The Angular Distribution of Sources 


Angular spatial methods attempt to determine if the observed 
distribution prefers a position in space (such as the Galactic Center) 
or a symmetry plane (such as the ecliptic plane, the Galactic plane, or 
the plane containing the Local Supercluster of galaxies) . Down to a 
minimum fluence of roughly 3 x 10"® erg/cm^, the spatial distribution 
of burst sources appears to be extremely isotropic (Mazets et al. 1981, 
Atteia et al. 1987). It also shows a negligible dipole moment, minimal 
quadrupole moments, and no propensity for clustering (Hartmann and 
Blumenthal 1988, Hartmann and Epstein 1989), in agreement with a random 
spatial distribution. 


2. The Radial Distribution of Sources 


a. Log (N) -Log (S) Curves 


The number of sources N brighter than some fluence S provides one 
method by which the radial (distance) distribution of the sample can 
be measured. The luminosity distribution of a sample of sources with 
spatial density n in the luminosity interval (L,L+dL) is described by 
the luminosity function <D(L) dL. Functionally, N(£S) - / n V <I>{L) 
dL, and, since the fluence of a source decreases as r^ (where r is the 
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distance to the source), S - L/(47ir 2 ). For a uniform spatial 
distribution of sources (where the volume V - [4rc/3] r ), integration 
yields N °c S“^^ r or log(N) « -3/2 log(S) . Similarly, a distance- 
limited sample will start off as log (N) «= -3/2 log(S), but will turn 
over at faint fluence such that log (N) « 0, A sample confined to the 
volume of a disk initially yields log(N) « -3/2 log(S); tilts over to 
become log(N) « - log(S) at fainter fluence, and finally reaches 
log (N) « 0 at faint fluence (when the disk has been completely 
sampled). A general description of log (N) -log (S) plots for disk 
models may be found in Fishman (1979) . 

Analysis of the log (N) -log (S) distribution of burst sources from early 
satellite catalogs (e.g. from KONUS, Mazets et al. 1981) indicates 
that the sample is in some way radially bounded, although this 
interpretation is quite controversial. 

Higdon and Lingenfelter (1986) have noted that KONUS looks for bursts 
on low-energy (soft) channels, whereas some bursts apparently have 
stronger high-energy (hard) emissions. They suggest that the observed 
sample is biased to brighter (nearer) sources (as both hard and soft 
sources are detected at high fluence, whereas hard sources are more 
difficult for the equipment to trigger on if they are farther away) , 
and argue that the corrected KONUS data is consistent with a 
log (N) -log (S) slope of -3/2. However, since the exact peak energy 
distribution (the relative number of soft bo hard sources) is unknown, 
the amount of correction applied to obtain these results is quite 
model -dependent . 

Since the detectors used trigger on minimum flux rather than on 
minimum fluence, they are more likely to sense short bursts than long 
ones . Although many observers suggest that the raw peak photon rate C 
should be used instead of fluence to bypass detector response, the 
overall problem of interpretation still remains unresolved. 

Paczynskii and Long (1986) infer that the faint-end log (N) -log (C) 
slope is -1.07 (indicating a radial limit in agreement with a Galactic 
disk model), while Jennings (1982,1984) corrects the data to a slope 
of -3/2. 

Meegan et al. (1985) have placed an upper limit on the log (N) -log (S) 
curve from balloon data. Their analysis incorporates (1) calculation 
of a detector count-rate triggering threshold, (2) conversion of this 
rate threshold to a fluence threshold, (3) simulations of the fraction 
of incident photons that will trigger each detector, (4) corrections 
for temporal triggering effects, and (5) analysis of the distribution 
based upon assumed spectral shapes of bursts. Their results indicate 
an upper limit of 2300 bursts/year at a fluence of 6 x 10" 7 erg cm" . 
In order to fall below this limit, the log (N) -log (S) curve must turn 
over below 10"^ erg cm"^ . 
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b. The Radial/Angular Distribution of Sources 


Schmidt (1968) developed another statistical test which analyzes the 
radial distributions of an astronomical sample, known as the V/V max 

test. This test is based upon the instrumental parameters C L £ m (the 
limiting count rate for a detector) and Cg (the source's observed 

count rate) . For a uniformly-distributed sample, the count rate of a 
source depends upon its distance R s . Since the limiting count rate 

corresponds the the distance Rmax that the source would have with 
count rate C Lim , K max « R <C s /C Lim ) 1 / 2 , and V/V^ - (C s /C Lim ) "3/2 . 

For a sample of sources, the average value <V/V max > « (Cg/C]^) “ 3 /2 

d (Cg/Cjjj^) “ 1/2. Thus, a sample with <V/V max > < 1/2 is distributed 

preferentially nearby, whereas one with <V/V max > > 1/2 is located 

preferentially far away. Schmidt, Higdon, and Hueter (1988) have used 
the V/V^x test on 13 HEA0-A4 gamma-ray bursts (with varying values of 

c Lim> to indicate that <V/V max > - 0.40 ± 0.08, which they point out is 

consistent with uniformity. However, this value is also consistent 
with the beginnings of a log (N) -log (C) turnover. 


C. Astrophysical Models of Burst Sources 


Many astrophysical models exist for gamma-ray burst production, and 
are summarized in reviews by Liang (1989) and Hurley (1989) . Possible 
models include compact sources in the cores of active galactic nuclei, 
massive galaxies, and globular clusters. However, the short timescale 
variability, cyclotron absorption features, and possible redshifted 
positron-electron annihilation lines in burst spectra lead to a 
favored model; that of a neutron star emitting by (1) accreting 
material and flashing, (2) vibration, (3) undergoing a crustquake, (4) 
being resurrected as a pulsar, or (5) via a magnetic flare. The 
limited power outputs of neutron star models suggest that a turnover 
in the log (N) -log (S) curve should be visible at some point. 
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Objectives 




Accurate models of the gamma-ray burst distribution are not obtained 
simply. Direct integrations of log (N) -log (S) and angular functions 
describing the distribution are generally quite difficult to perform. 
Also, the functional forms of these distributions are quite 
model-dependent, so direct integration of a great number of possible 
distributions is needed. This makes the problem even more unwieldy. 

Monte carlo techniques eliminate the difficulties of direct 
integration while simultaneously allowing .great freedom in examining 
model parameters. By creating a variety of discrete sources randomly, 
radial and angular distributions may be examined and compared to those 
of the actual data set. New models are easily created by merely 
adjusting model parameters. Models which are obviously incompatible 
with the observed data can be eliminated, while statistical tests 
performed on compatible models place restrictions on allowed parameter 
values . 

For this project, a Galactic disk spatial distribution is chosen for 
the sources (halo models can be easily incorporated in the future, 
although the log (N) -log (S) distributions of these models do not 
apparently turn over fast enough at low fluence) . A variety of source 
luminosity functions are tested within this framework, so that models 
with unreasonable luminosity functions might be eliminated. 
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A Monte Carlo Model of Galactic Disk Sources 


The procedure used to test the validity of a model is to (1) choose 
model parameters from which a monte carlo data set is generated, (2) 
select only those generated sources which have fluences bright enough 
to be measured, and (3) statistically compare the radial and angular 
properties of the artificial data set with the real one. A number of 
these models can be eliminated from consideration. Finally, general 
luminosity and spatial properties of the real data may be inferred 
from the monte carlo models. Creation of the monte carlo data set is 
carried out via program GAMMA. FOR (figure 1), and the data analysis is 
performed by program GAMMA. PRO. 


A. General Model Parameters 


The general model parameters are allowed to be changed within the 
program, and it is in fact not difficult to change the functional 
forms of any parameters. The gamma ray bursts are assumed to 
originate on Galactic neutron stars, so thd spatial distribution is 
that of the Galactic disk. Model parameters describing the Galactic 
disk are taken from the book Galactic Astronomy, by Mihalas and Binney 
(1981), and are meant to mimic general properties of the Milky Way (as 
presently understood) . These values may, however, be varied easily 
within the program. A variety of luminosity functions are also 
allowed in the model, and parameters describing these functions may be 
easily changed. 


1 . The Three-Dimensional Galactic Disk 


This disk model is intended to represent a variety of stellar ages and 
populations, but the youngest disk (spiral arm) population has been 
purposefully overlooked because the spatial distribution of bursts is 
apparently too isotropic to be associated with (1) spiral arms, and 
(2) a tendency to cluster as is observed for the youngest objects. 

The Galactic disk is therefore characterized 

(1) in the ©-direction by an isotropic distribution (i.e. there are 
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no angular structures such as spiral arms) . 

(2) in the z-direction (perpendicular to the Galactic plane) by a 
Gaussian distribution. There are physical reasons for choosing 
this model over either a linearly- or an 

exponentially-decreasing one (although both of these are often 
used in Galactic astronomy) : Primarily, these models exhibit 

discontinuities at z = 0, whereas a Gaussian model does not. 
Furthermore, these models are not as easy to explain in terms 
of maintenance via a Galactic gravitational potential. The 
width of the Gaussian distribution is represented by the 
variance CT Z 2 , where O z differs from the exponential scale 
height by ^2. Choosing a value for a z is difficult. Scale 
heights of stellar populations are estimated from Mihalas and 
Binney (1981) p. 278, to be as follows: Spiral Arm population 
(age 5 10 8 years) = 120 parsecs. Young disk population (age = 

10 8 years) = 200 parsecs. Intermediate disk population (age >=> 5 
x 10 9 years) = 400 parsecs, and Old Disk population (age ^ 10 10 
years) = 700 parsecs. A disk sample with a constant birthrate 
function would have an average age of around 5 x 10 9 years (if 
the disk is 10 10 years old) , and a scale height of roughly 400 
parsecs. This might be an underestimate, as the suspected 
scale height of white dwarf stars is 500 parsecs. A best guess 
value for o z is therefore 350 parsecs, corresponding to a scale 

height of 500 parsecs. 

(3) in the r-direction (radially outward from the Galactic center) 
by an exponentially-decreasing density function. This is both 
observed locally and from brightness distributions of other 
spiral galaxies of similar Hubble types. This exponential 
decrease is convolved with the infinitesimal radial increase (r 
x dr) needed to keep a constant disk density for cylindrical 
geometry. The exponential scaling factor Rg of the 

distribution exp(-r/R s ) is obtained from scaling the size of 
the Milky Way Galaxy to that of the Andromeda Galaxy (M31), and 
is therefore Rg =* 3.9 kpc . 

The sun is assumed to lie at a distance of 8.5 kpc from the Galactic 

center f in the Galactic plane. 


2 . The Luminosity Function of Burst Sources 


Several "standard” luminosity functions of burst sources are allowed: 

(1) All sources have the same luminosity, <L> . 

(2) The sources are chosen from a Gaussian luminosity function 
characterized by <L> and 
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(3) The sources are chosen from a topheavy luminosity function 

linearly increasing in number between and 

(4) The sources are chosen from a bottomheavy luminosity function 

linearly decreasing in number between and L^x* 

(5) The sources are selected from a bottomheavy power law 
luminosity function where N(L) « (L/I^^) 

A sample of the modeled disk distribution (for sources of a single 
luminosity <L>) is shown in figure 3. 


B . Selection of Data for Analysis 


Once each "burst” has been given a Galactic position and a luminosity, 
fluences and observed positions for it (in Galactic coordinates 1 and 
b) are calculated. Only those sources with fluences greater than S m ^ n 

are selected (S m £ n is chosen to be 2 X 10“ 7 erg cm"^ as this 

represents a minimum value of log(S) for which reliable log(N) data 
exists). From these data, a log (N) -log (S) array is built and angular 
characteristics are examined. An estimate must be made of the total 
number of sources, since this in part determines the "goodness of fit" 
of the log (N) -log (S) plot. At present this fit is estimated by 
inspection; future work may introduce a subroutine which optimizes the 
number of sources used. 

It should be noted that Galactic disk samples of low-luminosity 
sources take much more CPU time to create than those of 
high-luminosity sources. This is due to the sample's (disk 
height/limiting distance) ratio, which is larger for low-luminosity 
sources than it is for high-luminosity ones. In other words, the 
luminous sources are confined to an almost two-dimensional flat disk 
(which the computer rapidly fills) , whereas the low-luminosity ones 
are confined within a three-dimensional sphere (which fills in more 
CPU time) . 

Samples with a range of luminosities (such as the power law luminosity 
distribution described above) also use larger amounts of CPU time, 
which sometimes can be prohibitively large on a small computer such as 
a VAX . 
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C. Method of Data Analysis 


The log (N) -log (S) curve is produced from the data, and is chosen to 
(1) match the log (N) -log (S) plot observed from satellite data at high 
fluence, and (2) stay below the balloon limit (Meegan et al. 1985) of 
2300 bursts/year at a fluence of 6 x 10~ 7 erg cm -2 . The slope 
d [log (N) /log (S) ] /d [log (S) ] is also compared to that of the actual data 
set . 

The spatial distribution can be plotted for any minimum fluence in 
Galactic coordinates 1 and b. A quick-and-dirty analysis checks the 
distribution of events brighter than this minimum fluence as a 
function of Galactic latitude region by considering three latitude 
regions of equal area. The Galactic latitude regions 0° ^ |b| ^ 19.5° 
(low-latitude), 19.5° < \b I £ 41.8° (mid-latitude), and 41.8° < |b| £ 

90° (high-latitude), should contain equal numbers of stars n, with 
errors n/V(n-l). A more accurate method of analysis is also presented, 
as the system's dipole and quadrupole moments are calculated. A large 
dipole moment indicates that the distribution is biased towards one 
direction. The quadrupole moments yield two useful parameters: T| (the 
difference between the two closest eigenvalues of the quadrupole 
tensor) and £ (the most different eigenvalue) . When T| =* £ = 0 , the 
distribution is isotropic. An oblate (disk) distribution is indicated 
by £ > 0, whereas a prolate distribution (One which is strongly 
bipolar) is indicated by £ <> 0 . 


D. Results of Preliminary Data Analysis 


Samples of program output and data analysis are shown in figures 4, 5, 
and 6 . Once a number of models have been run, a comparison can be 
made between them and the actual data. 

The constraints imposed by both an isotropic angular distribution for 
high minimum fluence (S m j_ n £ 3 X 10”^ erg cm“ 2 ) and an apparent 

turnover in the log (N) -log (S) curve turn out to be quite strict. For 
o z ~ 400 parsecs and a single luminosity function L * <L>, a lower 

luminosity limit of <L> « 4 x 10 3 solar luminosities or 1 . 5 x 10 39 erg 
is necessary to turn the log (N) -log (S) curve over by S = 2 x 10~ 5 erg 
cm” 2 (the "elbow", or the lowest fluence where the log (N) -log ( S) curve 
can turn over in order to stay below the balloon data limit with a 
physically meaningful slope) . Models with lower average luminosities 
do not stay below the balloon data log(N) upper limit. Those with 
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higher average luminosities effectively stay below the balloon limit, 
but do not match the log (N) -log (S) curve at high fluence, do not show 
spatial isotropy for minimum fluences below that of the "elbow”, and 
predict a Galactic radius exceeding 15 kpc (compared to the apparent 
radius of 13.5 kpc). Thus only a small range of average luminosities 
between these extremes satisfies the observations, and the model is 
strongly constrained. 

Many implications arise from these limits: 

(1) The spatial density of sources is somewhere around 3 X 10“® 
pc”3 yr"l. If there are roughly 3 X 10^ neutron stars in the 
Galaxy (Hartmann, Epstein, and Woosley 1989), then the bursts 
sources must repeat on an average of around 10^ years to 
explain the burst rate. 

(2) The range of acceptable luminosities is too high for existing 
astrophysical models of neutron star crustquakes or resurrected 
pulsars (see Liang 1989), suggesting that these models are 
unlikely. 

(3) The effects of other luminosity functions on the results are 
not very pronounced. A Gaussian function broadens the "elbow" 
from a well-defined point to a curve and strengthens any 
spatial anisotropies that are present, due to an oversampling 
of the luminous sources (see also Hakkila 1989) . The topheavy 
and bottomheavy functions, being linear, exhibit slight effects 
that are very similar to the Gaussian distribution, although 
the bottomheavy function tends to show more spatial isotropy at 
high fluence. The power law function is the most promising 
luminosity function, but at present consumes an inordinate 
amount of CPU time. 
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Conclusions and Recommendations 


A. Conclusions 


This approach appears to provide a useful method for analyzing the 
properties of gamma-ray burst sources, as well as for analyzing 
properties of other astronomical objects. Detailed analysis which 
include interstellar absorption could lead to a better understanding 
of Galactic structure, stellar populations, and stellar evolution. 

Observations of gamma-ray bursts (angular isotropy at high fluence 
coupled with a log (N) -log (S) turnover) strongly constrain the number 
and types of allowed monte carlo simulations. The allowed models 
suggest that, for a distribution characterized by a single source 
luminosity <L>, only a small range of luminosities (10 39 erg ^ <L> 

< 6 x 10 ” erg) result in acceptable fits to the actual data. Other 

luminosity functions do not significantly alter these general 
conclusions, although the effects of bottomheavy power law functions 
are still unknown (due to prohibitive use of CPU time) . 

The implications of these results are worrisome, as they indicate that 

(1) some of the observational data is in error, (2) selection effects 
are still present in the existing data, (3) the Galactic disk model 
used in this monte carlo analysis is significantly in error, (4) a 
Galactic disk model is not the proper representation for the 
distribution of gamma-ray burst sources, or (5) the bursts really all 
have only a small luminosity range. 


B . Recommendations 


This approach is still only in the preliminary stages. More work 
still needs to be done in order to 

(1) quantify the comparison between models and observations, 

(2) vary disk parameters in order to see what effects different 
stellar populations have upon the resultant distributions, 

(3) examine the effects of a Galactic halo component, 

(4) integrate models over larger data samples to get better 
statistics, 

(5) include a method for optimizing the number of generated 
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sources, 

study the effects of a bottomheavy power-law luminosity 
function, 

(7) refine the Galactic model to incorporate new understanding 
about Galactic structure, 

(8) prepare for the possibility the BATSE will identify burst 
subpopulations, as these will then have to be examined 
separately, 

(9) etc. 

Of course, the new data obtained from the BATSE experiment on GRO 
should resolve the bulk of the burst distribution problem within only 
a few months of launch, and the data set that it generates can help 
isolate the locations of and mechanisms responsible for the gamma-ray 
burst sources. With the BATSE data, methods such as this will prove 
useful in isolating (1) the luminosity function of the burst sources, 
and (2) the stellar population of the burst sources (should they prove 
to have Galactic origins) . 
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Figure 1. The computer program GAMMA. FOR used to generate the 
monte carlo data sets. 
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2 • The computer program 
data sets. 
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ORIGINAL PAGE IS 
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mire 5 Sample output of monte carlo data analysis for the data 

Jet plottedTn t„. previous fig-- m ***** « 
the parameters used to generate the .model, the output 
latitude distribution of the sample, the dxpole «*»»«' 
parameters T\ and C of the quadruple moment for a V 

minimum fluences. The angular anisotropies present at ^ 
minimum fluence disappear when considering only bright source . 


XIII- 18 







